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Abstract 

We study a class of loop models, parameterized by a continuously varying loop fugacity n, on 
the hydrogen-peroxide lattice, which is a three-dimensional cubic lattice of coordination number 
3. For integer « > 0, these loop models provide graphical representations for «-vector models on 
the same lattice, while for « = they reduce to the self-avoiding walk problem. We use worm 
algorithms to perform Monte Carlo studies of the loop model for n - 0, 0.5, 1, 1.5, 2, 3, 4, 5 
and 10 and obtain the critical points and a number of critical exponents, including the thermal 
exponent y,, magnetic exponent y/,, and loop exponent y/. For integer «, the estimated values of 
yt and yi, are found to agree with existing estimates for the three-dimensional 0(n) universality 
class. The efficiency of the worm algorithms is reflected by the small value of the dynamic 
exponent z, determined from our analysis of the integrated autocorrelation times. 
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1. Introduction 

Among the many model systems studied in the field of statistical mechanics, two funda- 
mental examples that continue to play a central role are the n-vector model [1] and g'-state 
Potts model [2, 3, 4]. Initially, the parameters n and q can assume only positive integer values. 
However, the Kasteleyn-Fortuin mapping [5] transforms the Potts model into the random-cluster 
model [6], in which q appears as a continuous parameter. Likewise, certain «-vector spin mod- 
els on the honeycomb lattice can be mapped [7] to nonintersecting loop models which remain 
well-defined for non-integer n. Both types of mappings integrate out the original spin variables, 
while newly introduced bond variables define the remaining degrees of freedom. Each bond 
configuration can be represented by means of a graph covering a subset of the lattice edges. 
The transformed partition function specifies the statistical weight of each possible graph. This 
defines the probabilistic representation of the transformed model in terms of random geometric 
objects: clusters [5] and nonintersecting loops [7]. These geometric models play a major role in 



'Corresponding author 

Email addresses: liuqq@mail.ustc.edu.cn (Qingquan Liu), yjdengOustc . edu. cn (Youjin Deng), 
tim . garoniOmonash . edu (Timotliy M. Garoni), henkaiorentz . leidenuniv . nl (Henk W. J. Blote) 



Preprint submitted to Nuclear Physics B 



April 10, 2012 



recent developments of conformal field theory [8] via their connection with Schramm-Loewner 
evolution (SLE) [9, 10]. 

Given a particular lattice, or more generally a graph G = {V, E), we consider the loop model 
defined for n,x > Q by the partition function 

Z = ^«'-(^>xl'*l (1) 

A 

where the sum is over all configurations. A, of non-intersecting loops that can be drawn on the 
edges of G, \A\ denotes the number of occupied bonds, and c(A) denotes the number of loops. 
It is well known [7] that on any graph of maximum degree 3, the model (1) arises for positive 
integer n as a loop representation of an «-component spin model, 

Z = Tr +«xs, -S;), (2) 

ijeE 

where the s, = (s^, . . . , s") 6 R" are unit vectors and Tr denotes normalized integration with 
respect to any a priori measure (■)() on M" satisfying {s"s^}() = Safiln and {s")q — {s" sPs^)t^ - 0. 
In particular, uniform measure on the unit sphere is allowed, which results in an 0(«) symmetric 
spin model. In addition, various discrete cubic measures [7] are allowed, in which the spins 
are constrained to the unit hypercube; see Section 3.1. For each of these a priori measures, the 
model (2) with « = 1 is simply the Ising model. 

For n + \, the partition function (2) provides a high-temperature (small fi - xn) approxima- 
tion to the partition function of the n-vector model defined by the (reduced) Hamiltonian 

•H(Si,S2,...) = -y6^S,-S;. (3) 

ijeE 

For spins on the unit sphere, the Hamiltonian (3) with n - 2 and « = 3, defines the standard 
XY and Heisenberg models, respectively, while its n — » oo limit corresponds to the spherical 
model [II, I]. Since the partition function (2) shares the same symmetry as the Hamiltonian 
(3), one expects that for a given choice of spin a priori measure, the phase transitions of the 
two models will belong to the same universality class. We note that while (1) has a well-defined 
probabilistic meaning for all n, x > 0, the spin model (2) has positive weights only when x < l/n. 
Finally, we recall [12] that the « limit of the «- vector model corresponds to the self-avoiding 
walk (SAW) problem, however in this case the correspondence is not seen at the level of the 
partition function, but rather at the level of the two-point functions; see Sections 2. 1 and 3.1. 

On the honeycomb lattice, an exact analysis [13, 14, 15, 16] of the model (1) is possi- 
ble, which yields the critical point and the universal exponents as a function of « in the range 
-2 < n< 2. For three-dimensional lattices, however, comparable exact results are not available, 
and approximations are necessary. These exist in the form of renormalization at a fixed number 
of dimensions [17], the e-expansion [18], the l/n expansion [19, 20, 21, 22], series expansions 
[1] and Monte Carlo simulations [23]. See [24] and references therein. In particular, Wolff's 
embedding algorithm [25] is a highly-efficient Monte Carlo method for studying the 0(n) spin 
model for integer n > 0. 

Loop models with continuously varying n have recently been studied via Monte Carlo sim- 
ulations [26, 27] in two dimensions, however, we are unaware of any similar studies in three 
dimensions. In [27], worm algorithms were presented for simulating the loop model (1) on any 
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3-regular graph and for any real n > 0, and these algorithms were used to perform a systematic 
study of the loop model on the honeycomb lattice as a function of n. In the current article, we 
use the worm algorithms presented in [27] to perform an analogous study of the loop model 
on a three-dimensional 3-regular lattice, the so-called hydrogen peroxide lattice [28]. Since 
the hydrogen-peroxide lattice has coordination number 3, it provides a convenient setting for 
the study of high-temperature series [29, 30]. Therefore, given that worm algorithms simulate 
spaces of high-temperature graphs, it is also a natural setting for worm-algorithm studies in three 
dimensions. 

For our purposes, the hydrogen-peroxide lattice is perhaps best understood as a subgraph of 
the simple-cubic lattice. The vertex set is 7?, as for the simple-cubic lattice, however precisely 
half the edges are deleted from the neighborhood of each vertex, so that the edge set is Ej^UEyUE- 
where 

Ex = {{x,y,z){x H- z) : x,y, z e Z, z + x = Q mod 2), 
Ey = {{x, y, z){x, y + 1 , z) : X, y, z e Z, x + y = Q mod 2}, 
Ez = {{x,y,z){x,y,z + 1) : x,y,z e Z, y + z-0 mod 2}. 

Since this lattice is a cubic (i.e. 3-regular) subgraph of the simple-cubic lattice, a more descriptive 
name may be the "cubic-cubic" lattice. Being both 3-regular and transitive, it is in some sense the 
natural three-dimensional analogue of the honeycomb lattice. The minimal cycles have length 
10, with each vertex belonging to 10 such cycles and each edge to 15. A sketch of a finite 
patch wrapped on a 3-dimensional torus is shown in Fig. 1. Alternative drawings in which 
all angles are equal to 120° can be constructed [28] and are well-known in crystallography; the 
international number of such a lattice is 214 and its space group is /4i32. Both of these geometric 
configurations possess cubic symmetry. 

We shall employ two variants of the worm algorithm in our simulations. The first version 
directly utilizes worm transitions that are in detailed balance with the appropriate loop weights 
„e(A)^l>i|^ and is valid for any n > 0. The « — > limit of this version is equivalent to the 
Berretti-Sokal algorithm for SAWs. The second version combines the « = 1 worm algorithm 
with the coloring method [26, 27] to obtain an alternative algorithm which is valid for real « > 1 . 
The advantage of second algorithm is that it avoids the need to perform non-local connectivity 
queries. 

The details of these algorithms will be the subject of Section 2. Section 3 then describes 
the observables sampled in our simulations, which include a susceptibility-like quantity, and 
several dimensionless ratios. These are used to determine the critical points and the thermal and 
magnetic exponents by means of finite-size scaling. In addition, we present an estimator for 
the loop fractal dimension yi. Finally, Section 4 presents the results of our simulations and our 
conclusions are summarized in Section 5. 

2. Worm algorithms 

Given a finite graph G - (V, E), the cycle space, C(G), is the set of all A c E such that 
every site of G is incident to an even number of bonds in A. We call A c E and (V, A) Eulerian 



A /c-regular graph has precisely k edges incident to each vertex. The descriptions "/:-regular" and "coordination 
number i" are therefore synonymous. 
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Figure 1: A 64-site hydrogen peroxide (cubic-cubic) lattice with periodic boundary conditions. The sites occupy the 
same positions as in a 4 X 4 X 4 simple cubic lattice, however each site is adjacent to only 3 of the 6 edges allowed in the 
simple-cubic lattice. 



whenever A e C(G). The worm algorithm provides a natural method for simulating loop models 
of the form 

Z= ^ «^'-^>xl^l (4) 

AeC(G) 

where c(A) is the cyclomatic number of the spanning subgraph (V, A). The cyclomatic number of 
a graph is simply the minimum number of edges that need be deleted in order to make it cycle- 
free. Consequently, on graphs G of maximum degree 3 all Eulerian configurations A e C(G) 
consist of a collection of c(A) nonintersecting loops, and the models (1) and (4) coincide in this 
case. 

The essence of the worm idea is to enlarge the configuration space C(G) to include two 
defects (i.e. vertices of odd degree), and to then move these defects via random walk. When 
the two defects happen to occupy the same site they cancel, and the resulting configuration be- 
comes Eulerian once more. In the context of spin models, the state space of the worm dynamics 
corresponds to the space of high-temperature graphs of the two-point correlation function, as 
we describe in more detail in section 3.1. Worm algorithms for the high-temperature graphs of 
classical spin models were first formulated by Prokof'ev and Svistunov [31]. A careful investi- 
gation [32] of the dynamic behavior of the worm algorithm for the two- and three-dimensional 
Ising models showed that it suffers from only very weak critical slowing-down, especially in 
three dimensions where it was found that the Li-Sokal bound [33] appears to be very close to 
(or even exactly) sharp. In Section 4.5 we describe similar results for worm algorithms on the 
hydrogen peroxide lattice for n - 2. This contrasts markedly with the Swendsen-Wang algo- 
rithm, for which the Li-Sokal bound is far from sharp in three dimensions [34], and even more 
markedly with algorithms using simple local updates of the "plaquette" type [35], which have 
dynamic exponent z > 2. Furthermore, for certain observables, including the estimator for the 
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susceptibility, the worm algorithm was shown to display critical speeding-up [32]. A discussion 
of this estimator is the subject of Section 3.1. 

In [27], the worm methodology was used to construct algorithms for the loop model (4) on 
any graph G for any real « > 0, and these algorithms were used to perform a systematic study 
of the loop model on the honeycomb lattice. In the current article, we apply the algorithms 
presented in [27] to the hydrogen peroxide lattice. In this section, we give a brief description of 
these algorithms; refer to [27] for more details, including a discussion of a version which remains 
provably ergodic as x ^ oo. Related worm algorithms, corresponding to the loop expansions of 
both (2) and (3), are presented in [36, 37]. 

2.1. Simple worm algorithm 

It is most natural to describe the worm algorithm on an arbitrary finite graph G - {V,E). 
For simplicity, we assume G is regular (i.e. each site has the same number of neighbours). For 
distinct vertices «, v 6 V we let Su,v{G) denote the set of all subsets of E for which u and v are 
odd (i.e. have odd degree) with all other vertices being even, and we let .S,..,.(G) = C(G) for 
every v e V. The state space of the worm algorithm is the set S{G) of all ordered triples (A, u, v) 
such that M, V e V and A e S„^v{G). We emphasize that in the context of spin models Su,v{G) 
is nothing other than the set of high-temperature graphs of {s„ ■ s,,), so that S{G) corresponds 
to the space of high-temperature graphs of the magnetization; see Section 3.1. The proposed 
moves of the worm algorithm are very simple: starting from (A, m, v) e S{G) randomly choose 
one of the two defects; then randomly choose one of the neighbours of that defect; then move 
the chosen defect u to the chosen neighbour «' and flip the occupation status of the edge mm', 
so that A Aauu'. Here Aamm' denotes the symmetric difference of A with mm': i.e. delete 
mm' if it is occupied, or add it if it is vacant. One then simply applies a standard Metropolis [38] 
acceptance/rejection prescription to determine the acceptance probabilities required to ensure the 
worm transitions are in detailed balance with the desired configuration weights x''*' n''*'^' on S(G). 
The resulting Monte Carlo algorithm is summarized in Algorithm 1 . We emphasize that for any 
(A, u, v) € S{G), if M = V then A e C(G), so by only measuring observables when u - v one 
obtains a valid Markov-chain Monte Carlo method for the loop model. 

Algorithm 1 (Simple worm algorithm). 
loop 

With probability p choose the first defect, otherwise choose the second 
Uniformly at random, choose a neighbour u' of the chosen defect u 
Propose the update u ^ u' and A — > Aamm', leaving the unchosen defect v fixed 
Accept the proposed update with the acceptance probabilities given in Table 1 
if u' - V then 

Sample observables 
end if 
end loop 

The explicit acceptance probabilities used in our simulations are listed in Table 1 . Rather 
than the typical choice of min(l, z) for the acceptance function we chose z/(l +z), corresponding 
to a heat-bath algorithm [39]. We comment on this choice further below. We emphasize that, in 
general, the acceptance probabilities depend on the topology of the loops in a non-trivial way; 
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Table 1: Acceptance probabilities for the worm proposals (A, u, v) — ► (Aauu' , u' , v). The same acceptance probabilities 
are used for the proposals (A, i', u) — > (AAhu', i', u'). 



\Aauu'\ - \A\ 


c(Aamm') - c(A) 


P(Aamm',m',v)/P(A,m,v) 


Acceptance Probability 


+ 1 


+ 1 


nx 


nx/(l + nx) 


+ 1 





X 


x/(l + x) 


-1 


-1 


\l(nx) 


1/(1 


-1 





1/x 


1/(1 +x) 



we discuss this issue further below. The probability p in Algorithm 1 was set to 1 /2 in our 
simulations except when « = 0. 

The n — > limit of Algorithm 1 deserves some comment. By construction, the stationary 
distribution of the worm dynamics is 

h(A' ,11' y)eS{C) " 

As n 0, only those configurations with c(A) = will retain positive weight, and therefore 
the support of ncn.xi') in this limit reduces to the set of all paths (i.e. unrooted SAWs) on G. 
Likewise, the « limit of Algorithm 1 provides a valid Monte Carlo algorithm for simulating 
this reduced space of configurations. However, in the context of SAWs, one is generally inter- 
ested not in the set of all paths that can be drawn on G, but specifically in the subset of such 
paths which are rooted at a fixed site. By selecting p - \ m Algorithm 1 however, the location 
of the second defect becomes fixed and the — » 1, n — > limit of Algorithm 1 defines a valid 
Monte Carlo procedure for simulating the set of all rooted SAWs on G, with fugacity x for the 
number of steps. In fact, this algorithm is nothing other than (a slight variation of) the Berretti- 
Sokal algorithm [40] for the grand-canonical (i.e. variable-length) ensemble of SAWs. For the 
simulations reported in Section 4 therefore, when « = we used Algorithm 1 with p - I- 

An important practical matter when implementing Algorithm 1 is the need, when « + 0, 1, 
to perform a non-local query to determine if the cyclomatic number changes when an update 
is performed. Consider a spanning subgraph {V,A) c G. Since the number of components, 
k{A), is related to the cyclomatic number by k(A) - \V\ - |A| + c(A), the task of determining 
whether an edge-update changes the cyclomatic number is equivalent to determining whether it 
changes the number of connected components. The latter question can be answered by known 
dynamic connectivity-checking algorithms [41], which take polylogarithmic amortized time. A 
much simpler approach, which runs in polynomial time, but with a (known) small exponent is 
simultaneous breadth-first search [42]. We used the latter approach when implementing Algo- 
rithm 1 in the simulations presented in Section 4. In practice however it is worth noting that 
connectivity queries are not actually required at each step, as we now discuss. 

To illustrate, we will consider « > 1, and suppose the proposed update is to delete an edge 
mm' which is currently occupied. In this case the acceptance probability equals either 1/(1 -i- x) 
or 1/(1 + « x), depending on whether or not mm' is a bridge, and we have the bound 1/(1 + nx) < 
1/(1 + x) for all X > 0. In deciding whether or not to accept the proposed update, we draw a 
random number r uniformly from [0, 1]. If r < 1/(1 H- nx) then the proposal will be accepted, 
regardless of whether mm' is a bridge. Likewise, if r > 1/(1 + x) then the update will always 
be rejected. It is therefore only when 1/(1 H-nx) < r < 1/(1 -i-x) that we need to perform 
connectivity queries to determine whether or not mm' is a bridge, and therefore whether or not 
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to accept the proposal. A similar observation holds for the proposed addition of an edge, and 
analogous observations apply also for « < 1 . Finally, we note also that similar arguments can be 
used if the acceptance probabilities are chosen using the acceptance function min(l,z), however 
the results are rather more cumbersome in that case. This was our primary motivation for using 
the heat-bath version. 

Remark 2.1. It is possible to augment Algorithm 1 with the following additional move: when the 
two defects collide, move them both to a new vertex of G, chosen uniformly at random. It was 
found in [32] for the « = 1 case that the addition of this extra move does not alter the dynamic 
universality class of the algorithm, but does slightly improve its efficiency, with the efficiency 
gain tending to zero as the system size tends to infinity. For the simulations reported in Section 4 
we applied this extra move when n > 0. 

2.2. Coloring algorithm 

We now describe an alternative to Algorithm 1 which avoids altogether the need to perform 
connectivity queries when « > 1. Consider a finite graph G - {V,E). The key observation is 
that simulating the « > 1 loop model on G is equivalent to simulating the « = 1 loop model 
on appropriately chosen random subgraphs of G. To see this, let us begin with the loop model 
partition function (4) and introduce auxiliary vertex colorings as follows: 

Z= 2 n Z iScrc.r + (n-l)6^,.b], (6) 

AeC(C) Celcyclesof A| crce{r,h) 

" Z Zi n^-"-^ n ■^'^'[^<rc,r+(n-l)<5.c,b] n ^^..r. (7) 

AeC(G) [Te|r,b)'' ij<^ Ce|cycles of A) \'e(isolated vertices of A| 

Each vertex can take one of two colors, r(ed) or b(lue), and the leftmost product in (7) implies 
that each cycle has all of its vertices colored the same color. This fact allows us to unambiguously 
use the notation o-c for the color of a cycle C. The expression (7) defines a joint model of loops 
and vertex colorings, and one can simulate the loop model (4) by simulating the joint model (7). 
This is similar in spirit to the idea behind the Swendsen-Wang cluster algorithm [43]. Given a 
loop configuration A e C{G), the joint model (7) implies that we can obtain a random vertex 
coloring by coloring all isolated vertices red, and independently coloring each cycle red with 
probability l/n or blue with probability (n - !)/«. This coloring defines two random induced 
subgraphs of G, one red and one blue, and (7) implies that on the red subgraph the weight of 
each loop configuration is independent of «, and is in fact just the weight given by the n - I 
loop model. Therefore, the loop configuration on the red subgraph can be updated using an 
« = 1 worm algorithm. The weights of loop configurations on the blue subgraph do depend on «, 
but provided we re-choose the random vertex colorings every so often, to ensure ergodicity, we 
are free to simply leave the loop configurations on the blue subgraph fixed, since "do-nothing" 
transitions trivially preserve stationarity with respect to any measure. We therefore consider 
the red vertices as being active and the blue vertices inactive. These observations lead to the 
following algorithm for simulating the n > I loop model (4). 

Algorithm 2 (Colored worm algorithm). 
loop 

Current state is A e C{G) 
randomly with probability Pcoior do 
Color each isolated vertex red 
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Independently color each loop: red with probability 1/n, and blue otherwise 

Identify the induced subgraph of red vertices Gred = G[yred] 
end randomized block 
Choose, uniformly at random, v e V(Gred) 

Use n — I worm updates on Gred to make a transition (A^ed, v, v) — » (A^^^, v', v') 
New state is A' - A'^^^ U Abiue 
end loop 

The above arguments leading to Algorithm 2 can be expressed in a more formal setting using 
transition matrices, and a formal proof of validity can be given; see [27]. The underlying idea 
behind this coloring method is in fact very general [26]. A similar method was used in [44] to 
construct a cluster algorithm for the random-cluster model for arbitrary real ^ > 1 . Related ideas 
were also discussed in [45], in the context of a loop model with integer «. We also note that 
similar ideas were used in an exact mapping between an 0(«) and an 0(« - 1) model [46]. 

The probability pcoioi in Algorithm 2 can be set to any desired value in (0, 1]. We found 
empirically that in order to avoid a situation in which almost all of the computer-time is spent 
on re-coloring, rather than on worm updates, it is advantageous to choose pcoim to be strictly less 
than 1 . For each given value of n > 1 we therefore ran preliminary simulations to tune /Jcoior to a 
value which resulted in roughly equal time being spent of worm updates and color updates. 

2.3. Consistency checks 

We performed several tests to verify the correctness of the worm algorithms. First, we verified 
that Algorithms 1 and 2 produced the same results for a number of static observables, including 
the bond density. Furthermore, we checked the consistency of the results of the worm algorithms 
for general n with an alternative local algorithm which used plaquette updates. Some care is 
needed for these tests because, without special provisions, the plaquette algorithm is subject to 
conservation laws restricting the numbers of nontrivial loops spanning the periodic system. We 
also performed tests with cluster algorithms for the spin model (2), for the special cases « = 1 
and « = 2. For the latter case, the tests are restricted to the physical range x < 1/2, which just 
excludes the critical point, but still allows accurate tests. 

3. Observables 

In this section we define the observables that we sampled in our simulations, discuss some 
quantities of interest defined from these observables, and highlight some of their relevant prop- 
erties. The numerical results of our simulations and finite-size scaling fits for these quantities are 
then presented in Section 4. 

The following observables were sampled in our simulations. All observables were sampled 
only when the defects coincided, with the exception of the return time, T, which is defined on 
the full worm chain. The quantity L denotes the linear system size. 

• The number of bonds A/i(A) = \A\. 

• The number of loops Ni{A) - c{A). 

• The length of the largest loop Xi . 
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• The mean-square loop length 

where the sum is over all loops /. 

• The indicator "Wx for the event that there exists a loop with positive winding number in 
the X direction. We also sampled 'Wy and which are defined analogously. 

• The time T between consecutive visits to the Eulerian subspace of S{G) during a realiza- 
tion of the worm chain. 



From these observables we estimated the following quantities: 

• The expectations (Xi) and (Xa)- 

• The dimensionless ratio 

2 



The quantity Qi was used to locate the critical point Xc and estimate the thermal exponent 
y, when « > 0. It was not studied for « = since the configurations are single paths rather 
than multiple loops in this case. 



• The dimensionless ratio 

- (10) 

The quantity Qt, was used to locate the critical point Xc and estimate the thermal exponent 
jt when n = 0. It was not studied for n > since (A/^) is trivially of order L*'' for all x in 
that case. 

• The wrapping probability 

^(W.v + Wv+W.) (11) 

By symmetry, R - {'W^) - {"Wy) - (Wj), so the wrapping probability is simply the prob- 
ability of there existing a loop which has a positive winding number along a fixed coordi- 
nate axis. Analogous quantities were shown in [47] to provide a highly effective method for 
estimating the critical probability of square-lattice site percolation, and various probabili- 
ties of this kind can be calculated exactly for the random-cluster model on a square-lattice 
with toroidal boundary conditions [48, 49]. A notable feature of the wrapping probabilities 
studied in [47] was their weak finite-size corrections, and the results presented in Section 4 
show that such behavior also holds for the loop model on the hydrogen-peroxide lattice. 

• The expected time of return (7") to the Eulerian subspace. 

It can be shown for integer « > that (7~) coincides with the susceptibility x of the n- 
vector model (2), while for n = it coincides with the susceptibility of the self-avoiding 
walk problem; see Section 3.1 for details. The numerical results presented in [27] for the 
honeycomb lattice strongly suggest that in fact {T) ~ x holds near criticality for arbitrary 
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real n > 0^, regardless of whether one uses Algorithm 1 or Algorithm 2. Furthermore, the 
results presented in Section 4 suggest that (T) ~ x ^Iso holds on the hydrogen peroxide 
lattice for real non-negative n. 

We note that the quantities relating to loop properties, R, (-Ci), {££) and Q\, were studied for 
n > only. 

3.1. 'Worm return times and susceptibility 

In this section we provide a derivation of the result x - C^) for integer « > 0. We consider 
first the case n > 0. 

Let G = (V, £■) be a finite graph. The two-point correlation function of the model (2) is 



where the s, - (s^, . . ., s") e R" are unit vectors and Tr denotes the product over vertices / € V 
of normalized integration with respect to some given a priori (•)o measure on M". We will focus 
on three common cases for the spin a priori measure, corresponding to the 0(«), corner-cubic 
and face-cubic models. In each case, the a priori measure corresponds to uniform measure on 
a certain compact subset of W. In the case of the 0(«) model the spins are constrained to lie 
on the unit sphere; in the case of the corner-cubic model the spins point to the 2" corners of the 
unit hypercube; while in the face-cubic model the spins point to the midpoints of the 2« faces 
of the unit hypercube. We note that all three models coincide with the Ising model when n - \. 
Furthermore, when n - 2 the corner-cubic and face-cubic models coincide, since in this case 
they are simply related by a rotation of n/A. 

To identify [7] the spin partition function (2) with the loop partition function (4), one begins 
by expanding out the product over edges in (2), associating a bond of color a 6 {1, 2, ...,«} to 
each edge ij for which the term n x s" s" is selected. One can then attempt to explicitly sum 
over the spins. It is a common feature of each of the above a priori measures that the moments 
<(«')'"' (sj)"'- . . . (i")'"")o vanish unless each ma is even. Therefore, only terms for which each 
vertex is adjacent to an even number of edges of each color can make a non-zero contribution 
after performing the spin sums. One therefore arrives at a model of edge-disjoint colored loop 
coverings of G. For graphs of maximum degree 3, one can then simply sum over all ways of 
coloring the loops to obtain (4). 

The same general procedure can be applied to compute a graphical expansion of (12). In 
this case, the set of bond configurations leading to non-zero contributions is >S„,i,(G), the set 
of all bond configurations with precisely two odd vertices, u and v, introduced in Section 2.1. 
The summation over bond colorings of the component containing u, v (the defect cluster) must 
be treated somewhat carefully however, and it turns out that the precise result depends on the 
specific choice of a priori measure. In particular, different a priori measures can give different 
weights depending on the topology of the defect cluster. For graphs of maximum degree 3, the 
topology of the defect cluster is rather constrained and can be only one of the following: cycle, 
path, tadpole, dumbbell or theta graph. See Fig. 2. In general, whenever the defect cluster 




ijeE 



(12) 



^The case h = was not studied in [27]. 
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Figure 2: Possible topologies for the defect cluster. Paths and cycles are also possible, but not shown. 



is a theta graph it receives an additional model-dependent weight, in addition to x'"^' n'^^^\ A 
precise statement of the final result is given in the following lemma, whose detailed proof is 
straightforward and therefore omitted. 

Lemma 3.1. For any graph of maximum degree 3 we have 

Z(s„-s,)- 2 ;.l^l«^W0(A) 



6{A), if the component containing u and v is a theta graph, 
1 , otherwise. 




1 , face-cubic, 
(3«-2)/«~, corner-cubic, 
3/(« + 2), 0(n). 



We note that, unhke the corner-cubic and 0(n) cases, the result for the face-cubic model does 
not actually involve special weights for the defect cluster, because distinct colors correspond to 
orthogonal spin states in this case. Also, we note that, as expected, the expansion for the corner- 
cubic model coincides with that for the face-cubic model when n - 1,2. The 0(n) model only 
coincides with the face-cubic model for n - \. 

Now, from Lemma 3.1 we immediately obtain an expansion for the magnetization. By defi- 
nition M - Yiiev Si> so Lemma 3.1 implies 

<M'> = 2 ■ (^^^ 

= 1 2 xl^l«^(^'0(A), (14) 

{A.u.v)eS{G) 

where S(G) is the state space of the worm dynamics, as defined in Section 2.1, and Z is the 
partition function of the loop model (4). It should therefore come as no surprise that {A^) can 
be easily sampled via worm dynamics. 

Indeed, suppose we construct a worm algorithm on S(G) for which we chose the acceptance 
probabilities so as to give detailed balance with an arbitrary probability measure 



V^{A, u, v) oc ^{A), (A, u, v) € S{G). 
11 



(15) 



If denotes expectation with respect to P^(-), then it is obvious that 

<lc(C)V = 1^1 (16) 

where Zc(G)(iA) denotes the sum of over C(G), Zs(C){^) the sum of i^(-) over >S(G), and 1c(G)(') 
denotes the indicator for the event that (A, u, v) e S{G) satisfies A e C{G), or equivalently that 
u - V. Choosing (/'(A) = x^^^ rf^'^'' 0(A), and combining (14) and (16) while noting that states in 
C(G) do not contain theta graphs, it follows that 

(laa)> = ^ (17) 

where the expectation on the right is with respect to the «-vector spin measure defined by (2) 
for a given choice of a priori measure, and the expectation on the left is with respect to the 
corresponding graphical measure on the worm state space, as inferred from Lemma 3.1. The 
susceptibility of translation invariant systems-' is therefore equal to 

<Af)^^^ (18) 

^ |V| (lc(G,> 

In fact, the expression (18) also holds as a relationship between the « — > 0, p — > 1 limit of 
Algorithm 1, and the susceptibility of the self-avoiding walk problem. Indeed, taking this limit 
of Algorithm 1 implies 

<lc(c;)> = ^^ ^ ^. (19) 

ZjAe{rooted SAWs on C) ^ 

The denominator on the right-hand side of (19) is nothing other than the generating function 
of variable-length rooted SAWs on G, counted according to step length, which is precisely the 
definition of susceptibility for the SAW problem [50]. We note that in the « — > 0, /9 ^ 1 limit of 
Algorithm 1 the quantity \c(C) is simply the indicator for the 0-step walk. 

Finally, since (1c(G)) is the stationary probability of the worm chain occupying a state in 
C(G), the expected time between consecutive visits to C(G) is simply (7") - 1/(1c(G)); and so 
for integer n > we have 

{r)^x- (20) 

Remark 3.1. For the sake of computational efficiency, the simulations presented in Section 4 
used Algorithm 2 when « > 1 . Algorithm 1 and Algorithm 2 have precisely the same stationary 
distribution on the subspace C(G), and therefore the estimates they produce for all loop observ- 
ables defined on C(G) must be the same. However, they do give slightly different weights to the 
non-Eulerian states in S{G), which implies that (7~) will not be precisely equal to the x in this 
case. Nonetheless we expect that the scaling behavior of (7~) should still be governed by the 
magnetic exponent y/,. The results for the honeycomb lattice presented in [27] as well as our re- 
sults for the hydrogen peroxide lattice presented in Section 4 are consistent with this expectation. 
In order to test this expectation more systematically, we note that one can in principle construct 
a generalization of Algorithm 2 in which one introduces the vertex colorings to the full worm 
space <S(G) rather than just the subspace C(G). This would provide an interesting generalization 
of the coloring method presented in [26]. We shall investigate such methods elsewhere. 



'i.e. systems on regular lattices with periodic boundary conditions. 
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4. Simulations 



We simulated the loop model (4) on finite hydrogen-peroxide lattices with periodic boundary 
conditions, for loop fugacities n - 0,0.5,1,1.5,2,3,4,5 and 10. For n = 0,0.5,1 we used 
Algorithm 1, with p set to 1 for « = and to 1/2 for « = 0.5, 1. For « > 1 we used Algorithm 2, 
so as to avoid the need for non-local connectivity queries. 

For each value of «, we ran simulations of several LxLxL systems at several values of x. The 
length unit is chosen as one half of the 8-site hydrogen-peroxide cell, so that the systems contain 

sites. The lattice sizes L were taken in the range 8 to 128. Additional simulations for L - 256 
were performed at the critical point Xc as estimated from the smaller system sizes. For each 
choice of («, x, L) we performed multiple independent simulations, as a means of implementing 
(trivial) parallelization. For each independent simulation, statistical errors were estimated by par- 
titioning the data into 1000 bins, and computing error bars using blocking. The resulting means 
and error bars of the independent runs were then combined to produce the final estimates for each 
choice of («, x, L). For n + 0, we performed between 100 and 300 independent simulations for 
each choice of (n,x,L), each simulation consisting of 10^ returns to the Eulerian subspace, the 
first 10'* returns being discarded as burn-in. For n - 0, we performed 50 independent simulations 
for each choice of (jc, L), each simulation consisting of between 0.5 x 10^ and 1.25 x 10* 
iterations (hits), with the first 10^ L iterations being discarded as burn-in. For the n = case we 
sampled observables every /4 iterations, since we found in all cases that taking lags of this 
size guaranteed the autocorrelation function for Ni, was smaller than 0.1. (See Section 4.5 for 
relevant definitions concerning autocorrelations.) 

For each of the quantities discussed in Section 3 we performed a least-squares fit of our 
Monte Carlo data to an appropriate finite-size scaling (FSS) expression. As a precaution against 
excessive corrections to scaling, we imposed a lower cutoff L > Lmin on the data points admitted 
to the fit, and we studied systematically the effects on the fit due to variations of the value of Lmin- 
The error bars that we report for our fits are composed of the statistical error (one standard error) 
plus a subjective estimate of the systematic error. To estimate the systematic error we compared 
multiple distinct fits for each quantity, corresponding to different choices of correction terms 
included in the FSS formulas, and/or different choices of Lmin- Given that the FSS formulas are 
nonlinear, we used the Levenberg-Marquardt algorithm to perform the least-squares fits. 

4.1. The wrapping probability R 

According to finite-size scaling theory, the expected behavior of the dimensionless wrapping 
probability R{x, L) as a function of the linear system size L and temperature-like parameter x is 

R = + fli(x - Xc)U' + a2(x - x,fL^^'' + a^{x - x^^L^^'' + a4{x - x^fl'^^'' 

^ (21) 
+ c{x - XcfU' + biV + b2l^'- + biU' + e(x - x,)^'^^' + ... 

where x^ denotes the critical point, y, the thermal exponent, and yi ,y2,y3 are negative correction- 
to-scaling exponents. 

We fitted Eq. (21) to our Monte Carlo data, and the resulting estimates for x^ and y, are 
reported in Table 2. In each fit the values of the correction-to-scaling exponents were fixed, 
while Xc, yt were free parameters. For n - I, the correction-to-scaling exponents were fixed to 
yi = -0.815, y2 = -1.963, and = -3.375 [51], while for all other « > we set yi = -0.85, 
y2 - -1.8 and y^ = -3. The values of yi and y2 are based on estimates found in [52] for 
n = 0, 1, 2, 3 and the supposedly weak dependence of these exponents on n. 
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X X 

(a) R versus x for n = 1 . (b) Qi, versus x for n = 0. 

Figure 3: Wrapping probability R for the n = 1 loop model and dimensionless ratio Qi, for self-avoiding walks (n = 0), 
versus the bond (step-length) fugacity x, for linear system sizes 16 to 256. 



Table 2: Estimated critical point Xc and thermal exponent y, for a number of values of the loop fugacity n, as determined 
from least-squares fits of the wrapping probability R. 



n 


Xc 


.V/ 




a\ 


b, 


0.5 


0.5146078(6) 


1.651(5) 


0.1467(2) 


1.02(3) 


-0.08(3) 


1 


0.5180815(3) 


1.588(2) 


0.2524(3) 


1.41(1) 


-0.031(5) 


1.5 


0.5217878(5) 


1.539(4) 


0.3326(5) 


1.47(3) 


-0.062(5) 


2 


0.5257533(5) 


1.489(3) 


0.3990(4) 


1.48(1) 


-0.077(6) 


3 


0.5345904(8) 


1.398(3) 


0.5020(6) 


1.33(1) 


-0.144(5) 


4 


0.544904(2) 


1.329(8) 


0.5861(6) 


1.08(2) 


-0.18(2) 


5 


0.557096(5) 


1.275(12) 


0.6548(6) 


0.82(3) 


-0.21(2) 


10 


0.67367(2) 


1.142(15) 


0.8638(5) 


0.12(2) 


-0.23(2) 



For n < 3, each of the coefficients ai, . . . ,a^,bi,c were free parameters in the fits, while the 
coefficient e was set identically to zero. As a means of gauging systematic errors, we performed 
fits with b2, set identically to zero, as well as fits in which b2, b^ were free parameters, and 
compared the resulting estimates for y, and Xc- For each n < 3, our final estimates of y, and 
and their error bars were obtained by comparing these two fits. For « = 4, 5, 10, the strength of 
the corrections to scaling demanded that we include all terms in (21) in our fits, so that each of 
the coefficients oi , . . . , 04, , . . . , ^3, c, e were free parameters in these fits. 

4.2. The dimensionless ratios Qi and Qb 

As a consistency check on the values of x^ and y, determined from R, we performed anal- 
ogous fits of Qi, again on the basis of (21), with the same procedure for fixing the correction 
exponents and coefficient values. The resulting estimates for Xc and y, are reported in Table 3. 
The agreement with the corresponding estimates produced by the R fits is clearly excellent. 

The results for R and Qi applied only to n > 0, however analogous fits can be performed for 
« = using Qh- Following the same procedure as for 7? and Qi we thereby obtained the estimated 
values of x^ and y, for « = presented in Table 3. While this treatment of the SAW case appears 
natural within the context of the general loop model, we note that it may be considered somewhat 
unusual to deliberately study SAWs on a finite system. Indeed, one reason why SAWs are often 
considered easier to study via simulation than other models in statistical mechanics is precisely 
because it is possible to study SAWs on an infinite lattice. One might therefore be led to expect 
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Table 3: Estimated critical point and thermal exponent y, for a number of values of the loop fugacity n, as determined 
from least-squares fits of the dimensionless ratios Qi and Qt- 



n 


Xc 


yi 




a\ 







0.5113445(2) 


1.701(2) 


0.7017(2) 


0.675(4) 


0.11(2) 


0.5 


0.5146084(7) 


1.653(5) 


0.4858(6) 


0.578(6) 


0.012(6) 


1 


0.5180813(4) 


1.586(7) 


0.6158(4) 


0.836(5) 


0.042(5) 


1.5 


0.5217884(7) 


1.532(8) 


0.6884(5) 


0.750(6) 


0.052(5) 


2 


0.5257539(7) 


1.487(3) 


0.7316(4) 


0.644(6) 


0.066(4) 


3 


0.534595(4) 


1.398(4) 


0.7827(3) 


0.435(5) 


0.073(3) 


4 


0.544907(4) 


1.332(7) 


0.8086(4) 


0.270(6) 


0.080(7) 


5 


0.557108(8) 


1.280(10) 


0.8252(4) 


0.162(5) 


0.12(2) 


10 


0.67373(7) 


1.145(30) 


0.8573(3) 


0.016(5) 


0.10(2) 



Table 4: Estimated values of the scaling exponents y; and , for a number of values of the loop fugacity n, as determined 
from least-squares fits of (X2) and {T). For n < 1, Algorithm 1 was used, while Algorithm 2 was used for n> \. 



Fits for (£2) 


Fits for {T) 


n 


yi 




yii 


"0 


«l 









2.4875(7) 


0.796(5) 


-6.257(2) 


0.5 


1.723(3) 


0.64(2) 


2.482(4) 


1.11(2) 


6.8(1) 


1 


1.734(4) 


1.16(3) 


2.483(3) 


0.97(2) 


5.04(6) 


1.5 


1.755(3) 


1.44(3) 


2.482(3) 


0.91(2) 


3.80(5) 


2 


1.765(3) 


1.80(2) 


2.483(2) 


0.86(2) 


3.16(3) 


3 


1.795(3) 


2.18(2) 


2.482(2) 


0.80(2) 


2.32(4) 


4 


1.816(3) 


2.62(2) 


2.483(2) 


0.77(2) 


1.63(4) 


5 


1.834(6) 


2.9(1) 


2.483(3) 


0.74(2) 


1.22(4) 


10 


1.901(8) 


5.6(2) 


2.486(3) 


0.75(4) 


0.26(5) 



that Qh will be a rather poor estimator for the critical point. However, it is clear from Fig. 3(b) 
that this is not the case, and in fact the sharpness of the point of intersection of the Qt, curves 
in Fig. 3(b) is quite remarkable. It would be interesting to see if such behavior holds on other 
lattices, such as the simple-cubic lattice. 

4.3. Mean-square loop length 

The finite-size scaling behavior of {X.2) is expected to be [27] 

{£2) = L^^'^^ao + aiix - Xc)U' + biL>' + biU'- + b^U^ + ..) (22) 

where yi is the fractal dimension characterizing loop length. We fitted Eq. (22) to our Monte 
Carlo data for (-C2), fixing and Xc to their values as reported in Table 2. The corrections-to- 
scaling exponents were fixed according to the procedure described for R. And again, as a means 
of gauging systematic errors, for « < 3 we performed fits with b2, ^3 set identically to zero, as 
well as fits in which b2, b^ were free parameters, and compared the resulting estimates for y, and 
Xc- The resulting estimates for yi are reported in Table 4. 

4.4. Susceptibility X - average worm return time 

The finite-size scaling behavior of (7") is expected to be [27] 

<r) = l}^"'\aQ + ay{x - Xc)U'' + byU' + ..) (23) 
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where y/, is the magnetic exponent. We fitted Eq. (23) to our Monte Carlo data for (7"), fixing 
y, and Xc to their values as reported in Table 2 (Table 3 for n = 0). For the fits to (23) it was 
found sufficient to include only one correction-to-scaling term, and yi was fixed using the same 
procedure as for the R fits. The resulting estimates for yi, are reported in Table 4. 

A remark concerning our use of Algorithm (2) is in order. For « < 1, the values of (T) 
were computed using Algorithm 1, and so the discussion of Section 3.1 implies that for « = 0, 1 
the mean return time {T} is precisely equal to the susceptibility of the SAW and Ising models 
respectively. For n > 1 however, we used Algorithm 2 in our simulations, and so the discussion 
in Section 3.1 does not imply the exact equality of (T) andx in these cases. However, it was 
found empirically in [27] that on the honeycomb lattice while (T) may not be precisely equal to 
X when using Algorithm 2, the two quantities scale in the same way at criticality, so measuring 
{T} in Algorithm 2 was indeed found to be an accurate method for computing yi,. As we discuss 
further in Section 5, our estimates for yi, on the hydrogen peroxide lattice also agree within error 
bars with previously-obtained results for the 0(n) universality class. 

4.5. Dynamic behavior 

When n - I, Algorithms 1 and 2 coincide, and in this case a careful study of the dynamic 
critical behavior was carried out in [32]. In particular, it was found that in three dimensions the 
Li-Sokal bound (see below) appeared to be sharp. In order to explore the efficiency of the worm 
algorithms for more general n, we studied the dynamic critical behavior of both Algorithm 1 and 
Algorithm 2 for the case n - 2. For comparison, we also studied the dynamic critical behavior 
of a local plaquette update algorithm. 

Consider an observable O. A realization of the worm Markov chain gives rise to a time series 
0(t). The autocorrelation function of O is defined to be [39] 

{O(0)O(t)) - {Of 
Poit) = - 



where (■) denotes expectation with respect to the stationary distribution. From poiO we define 
the integrated autocorrelation time as 



1 °° 

Tint.O = 2 ^'5'^^^' ^'^'^^ 

and the exponential autocorrelation time as 

Texp o = lim sup ^— . (25) 

/^±co -log|po(OI 

Finally, the exponential autocorrelation time of the system is defined as 

Texp = sup Texp.O (26) 
O 

where the supremum is taken over all observables O. This autocorrelation time therefore samples 
the relaxation rate of the slowest mode of the system. All observables that are not orthogonal to 
this slowest mode satisfy Texp.o = T"exp- 

The autocorrelation times typically diverge as a critical point is approached, most often like 
T ~ where ^ is the spatial correlation length and z is a dynamic exponent. This phenomenon 
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is referred to as critical slowing-down [39]. In the case of a finite lattice at criticality, we define 
the dynamic critical exponents Zm.o^ Zexp.o and Zexp by 

Tint.O ~ -L-"'^ T^^p,0 ~ L''^"" Texp~L~°P. (27) 

Note that when defining Zexp and Zm.o it is natural to express time in units of sweeps of the 
lattice, i.e. U' hits. However, when using Algorithms 1 or 2, observables are sampled only when 
the chain visits the Eulerian subspace. Given that (T) scales like the susceptibility, this occurs 
roughly every L''^^^'' iterations, where Xt, is the magnetic scaling dimension. Since one sweep 
therefore takes of order L~-^'' visits to C(G), in units of "visits to CiG)" we have r ~ U'^^^i' . 
Consequently, we fit our dynamic data to the ansatz 

Ti„t,o = a + ^L-'-^+2x„ (28) 

to obtain the exponent Zint.o- 

In this work, we determined po(t) and Tint.o for the observables Nb, -Li, and A//. This was 
done for n - 2, using both Algorithm 1 and Algorithm 2. We find that the slowest of these 
observable is Nb- Figure 4 shows p^v,, as a function of f/Tint,//i, where a nearly pure exponential 
decay is observed, suggesting Zexp « Zim.TVfc- 
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Figure 4: Autocorrelation function PNf,(f) versus tlT\^f i^^^. The almost-pure exponential decay of PMi^(t) suggests 

For Algorithm 1 we find ZintM - 0-19 (3), while for Algorithm 2 we find ZmM - 0.243 (12). 
It is perhaps unsurprising that the bonds should decorrelate faster when using Algorithm 1 com- 
pared with Algorithm 2, since in the latter case a non-trivial subset of the bonds is always being 
held fixed at any given time. We emphasize, however, that these values of Zm.Ni, are inherent 
properties of the underlying Markov chains, and do not encode any information regarding the 
amount of CPU time required in practice to execute each MC update. In particular, despite the 
fact that Algorithm 1 suffers from slightly weaker critical slowing down, in practice Algorithm 2 
actually outperforms our implementation of Algorithm 1, because it avoids the need for costly 
non-local connectivity queries. 
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For comparison, we also determined Zint.Ni, for a simple plaquette-update algorithm. Like 
worm algorithms, these plaquette update algorithms can be used for arbitrary real « > 0. How- 
ever, a rough estimation for n - 2 yields Zmt,Ni, - 2-26 (5), which is significantly larger than 
for the worm algorithm. Similar values of z for plaquette algorithms have been observed previ- 
ously [35]. In addition, we note that the performance of the worm algorithm compares very favor- 
ably with Wolff's embedding algorithm [25], which is perhaps the most widely-used Monte Carlo 
method for studying 0(n) spin models, and which was found in [53] to have Zint.eneigy = 0.46(3) 
when (n, d) = (2, 3). 

Finally, it is interesting to compare these results with the n = 1 case [32] where it was found 
that Zint.yVfc ~ Q'/v ~ 0.17. This implies that the Li-Sokal bound [33] Zint.yv,, ^ o^/v^ which can be 
easily generalized to worm algorithms, appears to be sharp for « = 1 . By contrast, for n = 2 we 
have a/v x -0.026, so that the Li-Sokal bound does not appear to be close to sharp in this case. 

5. Discussion 

Using the worm algorithms introduced in [27], we have studied the critical behavior of the 
loop model (4) on a three-dimensional 3-regular lattice, the hydrogen peroxide lattice, for a range 
of integer and non-integer loop fugacities < « < 10. To our knowledge, this is the first direct 
Monte Carlo study of loop models with general n in three dimensions. A study of the dynamic 
critical behavior of these algorithms for n - 2 shows that the critical-slowing down is only very 
weak; not only is it much weaker than for local plaquette-update algorithms, but also significantly 
weaker than for the Wolff embedding algorithm [25, 53]. Combined with the previous study of 
the n = 1 case [32], these results strongly suggest that worm algorithms provide a very effective 
method for studying loop models in three dimensions. 

Our simulations show that the loop model (1) undergoes a continuous phase transition, with 
critical exponents that depend on the loop fugacity n. Our best estimates of the exponents y, and 
y/, are reported in Table 5. For integer n, these exponent estimates are consistent with the 0(n) 
universality class; for comparison, we list in Table 5 a number of relevant exponent estimates 
from the literature. For the cases « = and n - I this is entirely to be expected, since exact 
mappings unambiguously identify the loop model with the SAW and Ising models respectively, 
and these mappings are valid in a range of x which includes the critical point Xc- In fact, for the 
« = case Algorithm 1 actually directly simulates the grand canonical ensemble of SAWs. More 
generally, one might expect from the identity (2) that the loop model should in fact be in the 
same universality class as the n-vector model for any integer n. However, for n > 1 the situation 
is slightly more subtle, for two reasons. Firstly, the spin model to which the loop model maps 
has positive weights only for x < l/n, and for « > 2 we find empirically that Xc lies outside 
this range. Secondly, and perhaps more importantly, the mapping from the spin model (2) to the 
loop model (1) is many- valued; it maps the loop model (1) to a number of distinct spin models, 
whose Hamiltonians posses different symmetry groups and which therefore might be expected 
to belong to distinct universality classes. 

To put these observations in context, let us briefly recount the relationship between the 0(n) 
model and the two discrete cubic models discussed in Section 3.1. Renormalization group argu- 
ments [24] predict the following scenario. There exists a critical «£., believed to be ^ 3 in three 
dimensions, below which the cubic models and 0(«) models share the same universality class, 
but beyond which they do not. For n > n^, the face-cubic model is believed to undergo a first or- 
der transition, while the critical behavior of the corner-cubic model is believed to be governed by 
a distinct cubic fixed point. The discrepancies between the predicted exponent values of the 0(«) 
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Table 5: Comparison of the estimated loop-model exponents y, and y/, with existing literature. For h = we compare 
with known SAW exponents, while for integer n > 1 we compare with the corresponding exponents of the 0(n) spin 
model. No error bounds were given in [54] for their estimate of y, at n = 10. 



yi 


.V;. 


n 


This work 


Literature 


This work 


Literature 





1.701(2) 


1.7018(6)[55], 1.70179(5)[56], 1.70185(2)[57] 


2.4875(7) 


2.4849(5)[58] 


0.5 


1.653(5) 




2.482(4) 




1 


1.588(2) 


1.5868(3)[51] 


2.483(3) 


2.48180(8)[59],2.4816(1)[51] 


1.5 


1.538(4) 




2.482(3) 




2 


1.488(3) 


1.4888(2)[60] 


2.483(2) 


2.4810(1)[60] 


3 


1.398(2) 


1.406(10)[61], 1.405(2)[62] 


2.482(2) 


2.4830(5)[61], 2.481 1(15)[62] 


4 


1.332 (7) 


1.3375(15)[63], 1.333(4)[62] 


2.483(2) 


2.4820(2)[63], 2.4820(15)[62] 


5 


1.275(12) 


1.309(7)[64] 


2.483(3) 


2.4845(15)[64] 


10 


1.142(15) 


1.164(-)[54] 


2.486(3) 


2.488(1)[54] 



and cubic fixed points at « = 3 are too small to detect using current approximations, however they 
increase with n. Indeed, as n — » oo, the 0(n) model approaches the spherical model [11, 1], with 
exponents y/.spherkai - 1 and y/i^sphencai = 5/2, while the corner-cubic model can be reinterpreted 
as a constrained Ising model [65] so that its exponents are given by a Fisher renormalization of 
the Ising exponents, jf^cubic = ^/-JMsing = 1.4132(3) and, y/,,cubic = y/i.ising - 2.4816(1). For finite 
n these exponents have corrections of 0(1 /n); for example [66] predicts y/.cubic = 1.416(12), 
1.401(16), 1.404(12)for« = 3,4,8. 

Our estimates of y, reported in Table 5 agree within error bars with the corresponding 0(«) 
values from the literature when « < 4. For « = 5, 10 the agreement is not perfect, but still quite 
convincing. Furthermore, the «-dependence of our estimates is entirely consistent with the 
limiting 0(«) value of y, = 1 as n — > oo. The behavior of yi, in Table 5, while only weakly 
dependent on «, is also consistent with the 0(«) universality class. By contrast, our y, estimates 
are entirely inconsistent with the cubic fixed point. Based on these observations it seems reason- 
able to conclude that the critical behavior of the loop model ( 1) on the hydrogen peroxide lattice 
belongs to the 0(«) universality class. 

Finally, we note that the literature also contains results for the loop exponent yi in three di- 
mensions. Recently, Winter et al. [35] determined the fractal dimension of the high-temperature 
graphs at criticality for the Ising and the XY model on a cubic lattice. Using a plaquette update al- 
gorithm, they estimated the fractal dimensions of the Ising and the XY model as D - 1.7349 (65) 
and D - 1.7626 (66) respectively. Furthermore, Prokof'ev and Svistunov [67] reported the 
fractal dimension of the graph expansion of the complex \<p\'^ theory at its critical point as 
D - 1.7655 (20). From the expected equivalence in universality we may compare this results 
with our result of yi for the « = 2 model. As can be seen from the Table 4, our estimates of yi are 
in good agreement with the cited literature. 
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